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SUMMARY 


In this document, a method is developed for using the integrals of systems of 
nonlinear, ordinary, differential equations in a numerical integration process 
to (1) control the local errors in these integrals and (2) reduce the global 
errors of the solution. The method is general and can be applied to either 
scalar or vector integrals. A number of example problems, with accompanying 
numerical results, are used to verify the analysis and support the conjecture 
of global error reduction. 


INTRODUCTION 

Whittaker (ref. 1) defines the integral of a system of differential equations 
as a function of the state and time when the total time derivative is zero and 
the state variables are any functions of time. In solving nonlinear 
differential equations by a numerical integration process, these integrals will 
be used along with their associated numerical errors to (1) control the local 

errors in these integrals and (2) reduce the global error of the solution. 

I 

Integrals of systems of differential equations are constraints on the solution 
and, as such, can be used to reduce the number of degrees of freedom in the 
problem. Topologically, as more integrals of a system are introduced into the 
problem, the solution is constrained to lie on a larger manifold of the 
solution space with a resulting reduction in the global errors. The limiting 
case is, of course, an analytic solution to the differential equations where 
the solution is known at any time and the integral and global error are zero. 

The direct approach (analytical substitution) to using integrals of systems of 
differential equations to reduce the number of degrees of freedom in a problem 
and to reduce the errors introduces other types of difficulties (such as sin- 
gularities and switching logic in the remaining unsolved equations). 

Invariably, the overhead of calculating the right-hand side (RHS) of the 
remaining differential equations increases significantly; thus, any overall 
advantage of such an approach is nullified. 

In a direct analytical approach, Szebehely (ref. 2) linearized certain 
nonlinear differential equations by utilizing the integrals of the system and 
by introducing a new independent variable. This was an effort to formalize the 
results of Stiefel and Scheifele (ref. 3)> Burdet (ref. 4), Szebehely (ref. 5), 
and Sperling (ref. 6) who attempted to linearize, and in some cases stabilize, 
certain nonlinear differential equations by transformations of the dependent 
and independent variables of the problem. 

A direct numerical approach was made by Nacozy (ref. 7). The numerical errors 
in the integrals of the system were used to rectify the solution at each 
integration step in an attempt to stabilize the solution and reduce the errors. 
For some integrals, a linear expansion was necessary to compute the correction 
vector. Using a fourth-order predictor-corrector integration routine with a 
variable stepsize, global error reduction of two or three orders of magnitude 
was obtained. 



An optimization technique is the most general, indirect approach to controlling 
the error in a system described by nonlinear differential equations. A perfor- 
mance function is defined as a function of the error in the integrals, and the 
differential equations of the system are adjoined to this performance function 
by Lagrange multipliers. Through the use of variational calculus or some other 
similar technique, differential equations for the multipliers and optimality 
and boundary conditions can be developed. However, this is not an advisable 
procedure for controlling the errors for the following reasons: (1) the number 
of degrees of freedom in the problem typically increases twofold, (2) an it- 
eration procedure is introduced to satisfy the developed boundary conditions, 
and (3) the overhead in the solution of the problem increases significantly. 
However, the error in the integral would be minimized and would be indirectly 
used to reduce the global error of the solution. 

A functional, indirect approach to introducing integrals of systems in the solu- 
tion of differential equations was advanced by Baumgarte. In a number of 
studies (refs. 8 through 10) he used the integrals of the system to stabilize 
certain nonlinear differential equations and reduce the errors. The procedure 
was based on the principle of adjoining a form of the constraint (the coeffi- 
cient of the second derivative) by a Lagrange multiplier to the original sec- 
ond-order system of differential equations. The technique was applied to sev- 
eral problems and the results were encouraging. The error in the integral 
(almost always the energy) or constraint was substantially reduced by using 
this control process. However, several characteristics of these studies were 
disappointing: (1) global error results of the solution were almost always 

absent although some analytic solutions were available, (2) the technique re- 
quired a particular formulation for the problem, (3) the parameters that were 
introduced were not mathematically defined, and (4) a lack of generality 
existed when applying the approach to any system of equations and constraints. 

This document does not pretend to advance the best technique for using and 
controlling the errors in the integrals of a system of differential equations. 

It does propose a general technique for incorporating the integrals of a system 
of nonlinear differential equations and their associated errors in a process 
that will control the integral errors and reduce the global error of the solu- 
tion. The process requires no increase in the number of degrees of freedom in 
the solution. The technique has two disadvantages: (1) it requires some 

premathematical analysis to formulate the control vector, and (2) it generates 
some additional overhead and complexity in the solution. 


STATEMENT OF PROBLEM 

This study examines the numerical solution of a first-order, nonlinear system 
of ordinary differential equations of the form 


X = F(X,t) X(0) = Xq at t = 0 


( 1 ) 


which possess integrals of motion of the form 
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J(X,t) = K 


( 2 ) ■ 


where, X and F are n-vectors and J and K are m-vecbors; m < n and K 
is a constant defined by the initial conditions* The system of equation (1) is 
assumed not to be analytically integrable in terms of known functions and, hence, 
must be solved by a numerical integration process. Specifically, this method 
presents a numerical solution to equation (1) with the additional criteria 
that the solution lies "arbitrarily close" to the (m ■+• 1 ) -dimensional manifold 
described in equation (2). 

One obvious solution to this problem is to use equation (2) to eliminate m of 
the X*s and thus reduce the problem to integrating only the remaining (n-m)- 
first-order differential equations. However, past experiences and the examina- 
tion of a particular system of differential equations and their associated inte- 
grals indicate that this is not an advisable procedure. There is information 
about the solution embedded in equation (2), and it should not be used only as 
a check on the numerical integration of equation (1), but it should be used 
either directly or indirectly in the solution to reduce the errors and possibly 
the number of degrees of freedom in the problem. For example, if a problem- 
free procedure could be devised for introducing the m integrals of motion 
into a system of n differential equations, and the m integrals were intro- 
duced into the solution one at a time, the number of degrees of freedom in the 
solution would be reduced by one each time, requiring the solution to lie on 
a manifold with a dimension increasing by one. If there were actually n inte- 
grals of motion, then as m approaches n, the solution would be constrained ‘ 
to a larger subspace of the problem, and would also have global errors that 
tend to zero-vanish when m equals n. 

Because the direct approach to using the integrals of motion in the solution of 
a system of differential equations introduces other mathematical difficulties, 
an indirect approach shall be proposed to solve equation (1) while attempting 
to satisfy the constraints expressed by the integrals of the motion. 


SOLUTION WITH ERRORS 

If equation (1) is solved by a numerical integration process, the solution will 
certainly contain errors. Consider these errors as due to the inability of the 
numerical integration process to correctly evaluate the right-hand side (RHS) 
of the differential equations. Defining the numerical errors in the RHS ^ from 
the integration process as 6F, the differential equation (eq. (1)) can be 
expressed as 


X = F(X,t) + 6F 


(3) 


where at the initial time (t = 0), 6F(0) = 0. The vector F(X,t) in equa- 
tion (3) is the exact representation of the RHS of the differential equations. 
Now consider the term 6F as a perturbation to the original system of equation 
(1). The integrals of motion (eq. (2)) are not conserved (K is not equal 
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to a constant) and the differential equations for their rate of change can 
be expressed as 


K = G(X, 6F, t) 


( 4 ) 


where G is a vector function of the indicated arguments. The numerical error 
in the integral of motion is defined as 


£ = K - Ko Kq r K(0) 


(5) 


The differential equation for the time rate of, change of the error is simply 


e = K 


( 6 ) 


However, 6F is not known (except at the initial time) since the exact solution 
of the RHS in equation (1), F(X,t) is not known. Hence, it appears that this 
development is an interesting but insignificant exercise. 


SOLUTION WITH CONTROL 

In this section, a solution philosophy used in optimal control theory (ref. 11) 
will be adopted. In controlling a system, it is fundamental to have a process 
that is (1) observable and (2) controllable. The first criterion is certainly 
fullfilled, for it is noted in the numerical integration process that the value 
of the integral is not constant but grows in some manner characteristic to the 
particular numerical integrator, stepsize, etc. The second criterion, however, 
is not fullfilled. 

The process is not controllable because the error vector 6F is an unknown 
output of the numerical integration process and not an input. 

A control vector X (an n-vector) is added to the RHS of equation (1) (the 
exact equation) in an attempt to control the numerical error 5F of the inte- ■ 
gration process. 


X = F(X,t) + X 


(7) 


The justification for the addition of the vector X to the original equation 
(eq. (1)) is: if m integrals of the motion or constraints exist to a system 

of n equations, then the solution to this system of equations is constrained 
to lie on a m-dimensional manifold of the solution space and there are only 
n-m independent degrees of freedom in the problem. The introduction of a 
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control vector \ is an attempt to numerically reduce the n-order dependent 
system of equations to a (n-m)-order independent system of equations. 

The introduction of the control vector X to the RHS of the exact equation 
(eq. (1)) implies that the numerically integrated equation (eq. (3)) has a 
similar term added to its RHS. Of course, it is understood that with the 
addition of the control vector X, the term 6F appearing in equation (3) will 
have a different meaning since the error vector will certainly be disturbed by 
the introduction of this control vector. The introduction of the control 
vector X is an attempt to introduce a term to cancel .all or part of the error 
vector 6F such that when the numerical integration process is applied to equa 
tion (7), the error (defined in eq. (5)) will be arbitrarily close to zero. 
Thus, if 


6F + X = 0 


( 8 ) 


and equations (5) and (6) are used to obtain a "stable" solution for the 
control vector X, the numerical integration process will produce a value of 
the state X, which nulls the error. Also, since the number of degrees of 
freedom in the system of equations to be integrated have been numerically 
reduced, it is conjectured that the global error of the solution will also be 
reduced. A solution for the n-vector X from the system of m-constraint 
equations must now be obtained. 

Since the error rate and the error are now controllable, a stable differential 
equation for the desired functional relationship between these two errors is 
introduced as 


e = -Ye 


(9) 


where Y is a positive function (defined in the appendix). Now, any error 
arising in the integral of motion due to a dissatisfaction of equation (5) will 
be critically damped by the control vector X obtained from equation (9). 

Using equations (4), (5), (6), and (8) in equation (9) gives 


G(X,-X,t) = -Y (K(X,t) - Kq) 


( 10 ) 


Since the vector X must span-the-space defined by the state vector X, a 
solution for X is assumed to be 


X = AX 


( 11 ) 


where A is an undefined matrix of appropriate dimensions. Using equation 
(11) in equation (10) yields 
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G(X, -AX,t) = -Y (K(X,t) - Kq) 


( 12 ) 


The problem of controlling the integral errors is reduced to determining a solu- 
tion of an algebraic equation. Any solution for the undefined matrix A that 
satisfies equation (12) will produce a control vector (from eq . (11)) that, 
when used in the numerical integration process, will control the errors in the 
integrals of motion. The difficulty in determining a matrix A that satisfies 
equation (12) is, of course, problem dependent. In two of the example problems 
(linear oscillator and two-body problem) with a scalar integral of motion 
(energy), the matrix A was obtained by inspection. In a third example (two- 
body problem with an angular momentum integral), some manipulation was required 
to obtain a matrix A that satisfied equation (12). In a final example, a 
solution is developed to the two-body problem when both the energy and the 
angular momentum integral errors are present. 


ONE-DIMENSIONAL HARMONIC OSCILLATOR 

The first-order linear, differential equations describing the one-dimensional 
harmonic oscillator state are 


Xl = X2 

X2 = -X^ (13) 


Since an analytic solution to this problem exists, there are two integrals of 
motion. For this exercise, however, it is assumed that equation (13) cannot be 
solved analytically and that only one integral of motion (energy) exists and is 
defined as 

J(X) = 1/2 x^x = k 


where the superscript T refers to the transpose. 
The numerical error in the integral is defined as 


e = 1/2 X^x - ko 


ko = k(0) 


(14) 


Adding the control vector X to equation (13)» the controlled equation to be 
integrated is 


Xi = X2 + Xi 
X2 = -Xi + \2 
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( 15 ) 



Developing the total time derivative of the integral of motion, using equation 
(15) yields 


k = G(X,X) = X^X 


( 16 ) 


A stable differential equation for the functional relationship between the 
error and error rate is defined as 


e = -Y G (17) 

where Y is a positive scalar function. From equations (14), (16), and (17), 
the control vector X is required to satisfy the following equation. 


XT(X + - X) = Y ko 
2 ° 


( 18 ) 


If Y were a constant, equation (18) would represent a new integral of motion 
however, one that is functionally dependent on the control vector. 

A solution is assumed for the control of the form 



Yko aX 


(19) 


where a in an undefined coefficient. Using equation (19) in equation (18), 
a necessary condition is 


a X^x = 1 


One value of the coefficient satisfying the equation is 



which results in a control vector (from eq. (19)) of 
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X 


( 20 ) 


- - X 

2 k 


As expected, the control vector for the energy integral error is in the form of 
a feedback control law directly proportional to the error. It should be noted 
that the vector X is never singular unless the energy constant is zero 
(trivial case). 


Comparative Results 

Equation (15) was integrated with a fixed-step, fourth-order, Runge-Kutta 
integrator using the control vector defined by equation (20), Since the period 
of the uncontrolled solution is equal to 2tt, solutions were obtained for a con- 
stant integration stepsize defined as 


h = 2 tt/N 


where N is the number of steps in the integration process. The function Y 
was determined from a solution of the equation e(t + h) =0 at each 
integration step (see appendix). For this linear problem, the function Y for 
each value of N was found to be a constant. 

In the solution, it was noted that the uncontrolled error -in the energy grew in 
a linear manner with time. The controlled error remained essentially constant 
at a value five to ten orders of magnitude less (depending on the value of N 
and t) than the uncontrolled error. Thus, the method described in the analy- 
sis of controlling the error in the integral of motion appears to be valid for 
the linear problem. 

Because the analytic solution to' equation (13) is known, the global error of 
the numerically integrated solution can also be computed. The global error at 
a given time is defined as 

AX = |X^ - XA| (21) 


where the superscripts I and A on the vector X refer to the integrated 
and analytic values, respectively. The global error of the controlled solution 
was found to be always less than that of the uncontrolled solution. This indi- 
cates that the correct information from the error in the integral of motion was 
entering the solution via the control vector X. However, although the inte- 
gral errors were reduced substantially by the control, the global error showed 
only an infinitesimal reduction. This agrees with the results reported by 
Nacozy (ref. 7) for the harmonic oscillator problem. 
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Conclusion: For this linear and stable problem, controlling the error in the 

integral of motion has only a negligible effect on the global error of the 
solution. 


TWO-BODY PROBLEM WITH ENERGY INTEGRAL 

The first-order, nonlinear differential equations describing the two-body prob- 
lem are 


R = V (22a) 

V = -p ~ (22b) 

r3 


where R and V are the position and velocity vectors (respectively), ]s is 
the gravitational constant and r = | R| . There are three integrals of motion 
to this system of equations (not all independent) : two vector integrals - 
Laplace and angular momentum, and one scalar integral - energy. This example 
is concerned only with the energy integral that is formally obtained by scalar 

multiplying equation (22a) by V, and equation (22b) by R, then taking the dif- 
ference and noting the exact, differential* This integral can be expressed as 


J(R,V) = 1/2 vTv - p/r = k (23) 


where k is a constant of the motion defined by the initial conditions. If 
equation (22) is solved with a numerical integrator, the solution will contain 
errors. This numerical error in the energy integral is defined as 


€ = k - kQ k^ = k(0) 


(24) 


and the differential equation for its time rate of change is defined as 


• • 

e = k 


Now a control vector 

f Xy^) 
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is added to the RHS of equation (22) in an attempt to control the numerical 
error in the integral of motion* For convenience of notation, a state vector 
is defined as 



and an associated matrix is defined as 



where 0 and I are appropriate null and identity matrices, respectively. 
Using these identities, the numerical error and error rate can be expressed as 


£ = S^AS - ko 

£ = -X^S (25) 


The stable differential equation (eq~. (17)) is now used for the functional 
relationship between the error in the integral of motion and its time rate of 
change. Using equation (25) in equation (17) yields 


sT(X - Y AS) = -Yko (26) 

Since the vector X must span the space defined by the state vector S, an 
assumed solution for X is 


X - Y AS = -Yko A S 


(27) 


where A is an undefined matrix. Using equation (27) in equation (26), a 
necessary condition is 


S^AS = 1 


( 28 ) 
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Comparing equations (28) and (24), it is noted that one solution for the undefined 
matrix A is 



which results in a control vector (from eq. (27)) of 



(29) 


Thus, the control vector for the energy integral of the two-body problem can be 
cast in the same form as that of the control vector for the harmonic oscillator 
problem. It is only singular for the special case when the energy constant is 
zero (parabolic orbit). 


Numerical Results 

Equation (22) with the added control vector defined by equation (29) was in- 
tegrated with the same processor as the previous example. However, the con- 
stant stepsize was determined as 



where P is the period of the orbit defined by the initial conditions. 

Figure 1 illustrates a comparison of controlled and uncontrolled solutions for 
a two-body problem. The solutions were obtained using a fourth-order, Runge- 
Kutta integrator with a constant stepsize of N. The abscissa label of NORBIT 
refers to the number of orbits the equations were integrated, while the ordi- 
nate is the logarithm to the base 10 of either the absolute value of the posi- 
tion error or velocity error. (These two errors were approximately equal for 
this problem.) Hence, the ordinate is a measure of the number of significant 
digits of information remaining in the solution at a given time (a measure 
of the global error of the solution). The solid curves in figure 1 refer to 
solutions obtained with the addition of the control vector X, while the dashed 
curves refer to no control. 

In figure 1(a), the initial conditions were specified by a circular orbit that 
is an orbit with an eccentricity e = 0 and semimajor axis a = 1. For N = 20, 
the uncontrolled solution has lost all information content at NORBIT equals 20, 
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whereas the controlled solution still retains approximately two digits of ac- 
curacy. For N = 40, the controlled solution has approximately two additional 
digits of accuracy at NORBIT equals 40. For a given N, the two solutions are 
diverging; that is, the uncontrolled solutions have a steeper slope or are 
approaching zero more rapidly than the controlled solutions. 

In figure Kb) and 1(c), similar results are obtained for initial orbits with 
eccentricities of e = 0.1 and e = 0.2, respectively . For N = 20, the 
uncontrolled solution for e = 0.1 and e = 0.2 has zero-significant digits 
of accuracy at NORBIT = 15 and 9, respectively. 

When overlaying figures 1(a) and Kb), it is noted that the controlled solution 
shows almost no loss in accuracy for this change in eccentricity, whereas the 
uncontrolled solution shows significant losses. Similar (but not as 
significant) results are obtained when overlaying figures Kb) and Kc). 

The errors in the integral of motion (energy) in figures 1(a) through 1(c) 
behave similar .to that of the harmonic oscillator problem. The uncontrolled in 
tegral error grows almost linearly with time; the slope is determined by the 
initial conditions and the number of steps N. The controlled integral error . 
remains nearly constant in value, with a number of orders in magnitude less 
than the uncontrolled solution. The actual number of orders in magnitude 
depends on the initial conditions, integration stepsize, and integration time. 

In figure 2, the Y function is plotted versus time for one orbit for the 
three sets of initial conditions defined by e = 0.0, 0.1, and 0.2. For 
e = 0.0, which is similar to the previous linear problem, the Y function is 
nearly constant. The Y functions for e 0 are periodic in nature. For 
e = 0.2, there is a portion of the solution where Y < 0. Since there are no 
constraints imposed on either the sign or the magnitude of the function Y, 
a negative value is certainly possible (although disconcerting) when viewed 
from a solution to equation (9). 

In the next two example problems (concerned in part with controlling the error 
in the angular momentum vector), only the control vector and an outline of the 
solution are developed. 


TWO-BODY PROBLEM WITH ANGULAR MOMENTUM INTEGRAL 

The differential equations to be integrated are the same as in the previous 
example; however, the angular momentum integral shall now be considered, and it 
is expressed as 


J (R,V) = R X V = K (30) 


where the vector K is a constant of the motion defined by the initial con- 
ditions. Similar to the previous example, the numerical error vector in this 
integral is defined as 
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Kq = K(0) 


(31) 


e 


K - Kr 


T T 

and after introducing a control vector = (Xp , Xy) into the RHS of 
equation (22), the time rate of change of the error may be determined as 


e = R X Xy - V X Xp 


(32) 


A stable vector differential equation for the functional relationship between 
the vector error and its rate is defined as 


G = - “Y e 


(33) 


where Y is again a positive function. The control vector X must now span 
a space defined by the vectors R, V, and K or 


(34) 

where the ot»s and 3 's are undetermined scalars. But from equation (32), 
any component of the vector Xp parallel to the vector V and, similarly, any 
component of the vector Xy parallel to the vector R will have no 

effect on the error rate vector g. 

Hence, the vectors Xp and Xy have the simpler form 
Xr = ap R + K 

Xy = Bv V + K (35) 

Using equation (35) in equation (32), the error rate vector is 


X = 

OlR 

ay 

a^ 


R 

V 



^R 

1 

QQ. 


K 


e = (Or + By) K + (Bk R - Or V) X K 


(36) 


The first and second terms on the right-hand side of this equation are the rates 
of change of the angular momentum error vector along and perpendicular to the 
vector K, respectively. Taking the inner product of the vectors K, U = R x K 
and W = V x.K (respectively) with equation (36) and using equation (33) yields 
the necessary conditions on the coefficients a and B : 
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( 37 ) 


(apj + B\j) - -Yq 

Bk U^U - OLfr U^W = -Y uTe 

Bk U^w - W^w = -Y w^e 


where 


eT K 

kT~k 


The first equation implies that ap and By are homogenous in the term q. 
One solution for these coefficients is 


Op = -Y c q , By = -Y(1 - c)q 


(38) 


where the coefficient c is to be determined. The latter two equations (37) 
can be solved for and B^ as 


«K 


- KeTwHU^u) - (eTu)(uTw) 

A L 


3k 


= ^ j^(eTw)(uTw) - (eTu)(wTw) 


(39) 


where the determinant, A, is 
A 5 (uTu)(W^W) - (uTw)2 

It should be noted that the determinant is only zero for rectilinear motion. 

The control vector for the vector angular momentum error is 

Xp s cxp R + oijj K 

(40) 

Xy = By V + Bp K 
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where a and B are defined by equations ( 38 ) and (39) and Y is determined 
from the condition je(t + h)| =0 (see appendix). 

The coefficient c in the angular momentum error control vector is a scaling 
parameter between the components of error along and normal to the vector K 
(eq. ( 36 )). Since the physics of the problem dictate that the error normal to 
the vector K will be extremely small, the actual value of the coefficient c 
used in the solution will have a minimal effect on the global error. One 
difficulty should be mentioned. Since the angular momentum error is formed by 
a vector product rather than as in the energy error (a scalar product) a con- 
siderable amount of algebraic manipulation will be required to obtain the 
function Y. 


TWO-BODY PROBLEM WITH ANGULAR MOMENTUM AND ENERGY INTEGRALS 

In this section the scalar energy integral is added to the vector angular momen- 
tum integral, and a control vector X for both of these errors is determined. 
From equations (28) and (37), the necessary conditions on the coefficients a 
and 3 are 


yo^R 

r 


+ 6v VTV = 


aR + 6v = '^K 


(41) 


where aj( and 3 k given by equation (39) and the subscripts k, K refer 
to variables associated with the energy and angular momentum integrals, respec- 
tively. Solving equation (4l) for the coefficients aR and 3v yields 


OR = \ vTv)/(2k + y/r) 

e ll^ 0 


(42a) 


Bv = \ l^/r)/(2k + y/r) 


However, this solution is singular when the eccentricity of the orbit e = 0. 
A solution valid for e = 0 is 


« 


OR 


- Yk £k 
k 


6v 


Yk £k 
2k 


e = 0 


(42b) 
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with the additional condition 


2k 


= y 


K ^ 


( 43 ) 


Thus, for orbits with e = 0, the function Tp- is not directly determined by 
the integrator but by the function and the integral errors. But, the coef~ 

ficients given by equation (42b) will produce the energy error control vector 
given by equation (29). Then if the coefficients and are small and 

if e 0, the control vector that was used to reduce the energy integral error 
will also reduce the angular momentum integral error. The conclusion has been 
numerically verified. 

For the problem illustrated in figure 1(a) with N = 20, the energy integral 
control vector reduced the angular momentum integral errors by approximately 
five orders of magnitude. However, when the eccentricity e was equal to 0.1, 
the energy integral control vector only reduced the angular momentum integral 
errors by approximately one order of magnitude. 


CONCLUDING REMARKS 

This method for using the integrals of systems of nonlinear differential equa- 
tion and their associated numerical errors should be applied to (1) a variable 
step integrator, preferably the Runge-Kutta 4/5, (2) other problems where inte- 
grals or other type of constraints are satisfied through rectification at each 
integration step, and (3) an unstable system of nonlinear differential 
equations. 


Lyndon B. Johnson Space Center 

National Aeronautics and Space Administration 
Houston, Texas, March 28, 1980 
910-44-00-00-72 
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Figure 1.- Continued. 
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APPENDIX 


CALCULATION OF THE FUNCTION Y 


For a fourth-order, fixed-step Runge-Kutta integrator, the state vector is 
updated at a time t + h by the expression 


X(t + h) r Xq + - (f1 + 2(f2 + f3) + f'*) 
6 


Xo = X(t) (A1) 


where h is the integration step and the functions fJ are defined as' 


= F(Xo, to) 

f2 = f(xo ^ I Fl. to . 

f3 = f(xo + ^ f 2, to f ) 
= F(Xo + hF3, to + h) 


For the controlled solution, the right-hand sides of the differential equations 
may be separated into two parts; the vector F and the vector n adjoined by 
the unknown function 

X = F(x,t) + y n(X,t) 


where Y q = X. For simplicity of notation, the argument t is excluded from 
the vectors F and q. Then the vector may be expressed as 

Fl = F(Xo) + Yq(Xo) = Fq + Y qo (A2) 

and the vector f 2 may be expressed as 


A-1 



P'2 = F [Xo + + Y no)) + Y n (Xo + -(Fo + Y Ho)) 

2 2 


But the vector X represents a small perturbation to the vector F, and 
hence the vector may be approximated as 

p2 = Fi + Y (n1 + 6i) + ... (A3) 


where 


Fi = F(Xi) ni = n(Xi) Xi = Xo -H 5 Fo 


and 


<Sl 


^ (^ + y 

2 \3X ^ 9X/ Xi 


6o = 0 


Similarly, the vectors f 3 and F^ may be expressed as 

f3 = F 2 + Y(n2 + 62) 

F^ = F3 + Y(n3 + 63X 


where 


X2 = Xo + ^ Fi 

X3 = Xo + h F2 


Thus, the coefficient of 


in equation (A1) may be expressed as 


p) + 2(p2 + p3) + F^ = *? + Y(^7 + *f ) 


iM) 


(A5) 


A-2 



where 


= Fq + 2(F^ + F2) + F3 

"57 = Tio + 2(n1 + 7 i 2) + n3 

^ = 60 + 2(6'| + 62) + 63 

For clarity, the analysis is restricted to a particular problem - the harmonic 
oscillator. Extensions to other problems are straightforward. The desire 
is to determine the function Y, which will result at a time t + h in an 
integral error of zero. 

0 = e(t + h) = 1/2 X(t + h)T X(t + h) - ko (A6) 

Using equations (Al) and (A5) in equation (A6) yields ' 


Xo + g [/+ + 


= 2ko 


(A7) 


This is a polynominal in the function Y; however, if the small term — 

3X 

is ignored, equation (A?) may be expressed as the quadratic equation 
ay2 + bY+c = 0 

where 



0 = I X + 5 1^ - 2ko 

0 


A-3 



There are two solutions for the function however, it may be readily 
verified that the plus sign of the radical is correct. 
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